SOAP Séance 1: analyse de Fourier, échantillonnage

Table of Contents

• 1  Premiers pas: notion de signal et exemples
  • 1.1  Notion de signal
  • 1.2  Les sons en pratique
    • 1.2.1  Jeux de données
    • 1.2.2  Charger et afficher un son
    • 1.2.3  Fréquence d'échantillonnage
    • 1.2.4  Sauvegarder un son dans un fichier
  • 1.3  Les images en pratique
    • 1.3.1  Jeux de données
    • 1.3.2  Charger, afficher, sauvegarder une image
  • 1.4  Conclusion de la partie
• 2  Introduction à la transformée de Fourier discrète
• 3  Notion de fréquence
  • 3.1  Signaux périodiques
  • 3.2  Sinusoïdes
  • 3.3  Les sinusoïdes complexes
  • 3.4  Notion de fréquence normalisée/réduite
  • 3.5  Conclusion de la partie
• 4  Analyse spectrale: la transformée de Fourier discrète
  • 4.1  Représentation fréquentielle: intuition, exemples
    • 4.1.1  Représentation fréquentielle d'une sinusoide
    • 4.1.2  Représentation fréquentielle d'un son
    • 4.1.3  Les signaux périodiques et leur spectre harmonique
  • 4.2  Transformée de Fourier discrète (DFT)
    • 4.2.1  Intro: différentes transformées pour différents cas
    • 4.2.2  Définition et propriétés de la DFT
    • 4.2.3  Illustrations numériques des propriétés de la DFT
      • 4.2.3.1  Base orthogonale
      • 4.2.3.2  Périodicité et symétrie
    • 4.2.4  Représentation alternative avec fréquences négatives
    • 4.2.5  Algorithme rapide: fast Fourier transform (FFT)
  • 4.3  Conclusion de la partie
• 5  Fréquences spatiales et DFT pour les images
  • 5.1  DFT 2D d'une image
  • 5.2  Spectre 2D d'une image
  • 5.3  Spectre 2D d'un extrait d'image
  • 5.4  Les fréquences horizontales
  • 5.5  Les fréquences verticales
  • 5.6  Les fréquences obliques
  • 5.7  Conclusion de la partie
• 6  Comment échantillonner un signal ou une image?
  • 6.1  Les fréquences discrètes dans la DFT
  • 6.2  Théorème d'échantillonnage de Shannon-Nyquist
  • 6.3  Formule de reconstruction
  • 6.4  Comment sous-échantillonner

# Pour que les changements dans les modules importés (par exemple soap_utils.py) soient pris en compte
%load_ext autoreload
%autoreload 2
# Pour afficher les figures dans le notebook
%matplotlib inline

import warnings
from matplotlib import MatplotlibDeprecationWarning
warnings.filterwarnings("ignore", category=MatplotlibDeprecationWarning)

Utiliser le code de la cellule ci-dessous pour désactiver le scrolling quand l'affichage est long.

%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

import numpy as np
import matplotlib.pyplot as plt

import matplotlib as mpl
mpl.rcParams['figure.figsize'] = 20, 4
mpl.rcParams['lines.linewidth'] = 3
mpl.rcParams['font.size'] = 20
mpl.rcParams['lines.markersize'] = 15
mpl.rcParams['axes.grid'] = True
# mpl.rcParams.find_all('')

Premiers pas: notion de signal et exemples

• que signifie échantillonner? qu'est-ce qu'un échantillon?
• qu'est-ce qu'une onde?

Notion de signal

Un signal est un vecteur $\mathbf{x}\in\mathbb{R}^N$ représentant une séquence de coefficients appelés échantillons. $\mathbf{x}[n]$ est l'amplitude du $n$-ième échantillon. Exemple: un son est un signal.

La fréquence d'échantillonnage (sampling frequency) est le nombre d'échantillons par seconde. Exemple: $f_s=44100$Hz pour un CD.

Questions:

• quelle est la durée d'un signal de taille $N$ et de fréquence d'échantillonnage $f_s$?
• quelle est l'expression de l'instant $t_n$ du $n$-ième échantillon en fonction de $n$ et de $f_s$?

Les sons en pratique

Jeux de données

Jeu de sons de référence: SQAM

• Format .wav, stereo
• 70 sons (ajouts personnels possibles!), contenus variés (voir fichier pdf descriptif)
• Disponible sur Ametice, à installer dans data/sons/

Jeu de données du projet:

• sons de parole, plusieurs locuteurs, 1 fichier audio par locuteur
• disponibilité à venir

Charger et afficher un son

from scipy.io import wavfile  # Pour lire et écrire des fichiers sons

sound_path = '../../data/sons_int/'
sound_filename = sound_path + '35.wav'
fs, x = wavfile.read(sound_filename)
print(f'Son chargé: type {type(x)}, dimensions {x.shape}, dtype {x.dtype}')
print(f'fs = {fs} Hz')
print(f'Durée du son: {x.shape[0] / fs} secondes.')

Son chargé: type <class 'numpy.ndarray'>, dimensions (2601900, 2), dtype int16
fs = 44100 Hz
Durée du son: 59.0 secondes.

from soap_utils import plot_sound  # Pour afficher un son

plot_sound(x, fs)
plt.legend(['Left', 'Right'])

<matplotlib.legend.Legend at 0x7f7c3b1895e0>

Fréquence d'échantillonnage

(rappel) Fréquence d'échantillonnage fs = nombre d'échantillons par seconde.

print(f"Fréquence d'échantillonnage: {fs} Hz")
print(f"Période d'échantillonnage: {1/fs} seconde.")

t0 = np.round(0.5715 * fs) / fs  # Start time for zoom
plot_sound(x, fs, marker='.')
plt.xlim([t0, t0+1e-3])  # Zoom
plt.plot([t0 + 1 / fs] * 2, plt.ylim(), ':')
plt.plot([t0 + 2 / fs] * 2, plt.ylim(), ':')

Fréquence d'échantillonnage: 44100 Hz
Période d'échantillonnage: 2.2675736961451248e-05 seconde.

[<matplotlib.lines.Line2D at 0x7f7c14755970>]

Sauvegarder un son dans un fichier

y = x[:15*fs, 0]
plot_sound(y, fs)
wavfile.write(filename='excerpt.wav', data=y, rate=fs)

Astuce. Attention à l'encodage des sons! En général, les coefficients sont soit de type int16, soit float. Dans ce dernier cas, ils sont normalisés, entre -1 et 1. Mais si vous convertissez des coefficients de type entier en flottant et que vous sauvegarder un son sans le normaliser (c-à-d qu'il y a des flottants plus grand que 1 en valeur absolue), cela créera des fichiers sons distordus qui risquent de faire mal aux oreilles (et au matériel) lors de l'écoute. Recommandation: avant de sauvegarder un son, assurez-vous qu'il est de dtype entier ou que ses valeurs sont entre -1 et 1.

Les images en pratique

• Image en niveau de gris de hauteur $H$ et largeur $W$ $\equiv$ une matrice $\mathbf{M}\in\mathbb{R}^{H \times W}$.
• Élément $\mathbf{M}[l, c]$ $\equiv$ niveau de gris du pixel à la ligne $l$ et la colonne $c$.
• Le noir est codé par $0$; le blanc est codé par une valeur maximale, par exemple 1 (coefficients flottants) ou 256 (coefficients entiers)
• Image couleur de hauteur $H$ et largeur $W$ $\equiv$ tenseur $\mathbf{M}\in\mathbb{R}^{H \times W \times C}$ où $C$ est le nombre de canaux.
• Plusieurs codages de couleur existent, sur 3 ou 4 canaux. Par exemple Rouge/Vert/Bleu (RGB) sur 3 canaux.

Jeux de données

Jeu d'images de référence

• Format .png, niveaux de gris
• 19 images (ajouts personnels possibles!), contenus variés
• Disponible sur Ametice, à installer dans data/images/

Vidéos du projet

• disponibilité à venir
• possibilité d'extraire les images des vidéos

Charger, afficher, sauvegarder une image

image_path = '../../data/images/'
image_filename = image_path + 'baboon.png'
img = plt.imread(image_filename)
print(f'Image chargée: type {type(img)}, format {img.shape}, dtype {img.dtype}')

Image chargée: type <class 'numpy.ndarray'>, format (512, 512), dtype float32

plt.imshow(img, cmap='gray')  # Affichage

<matplotlib.image.AxesImage at 0x7f7c07901a60>

img2 = img[:100, 100:250]
plt.imsave(fname='eye.png', arr=img2)  # Sauvegarde
plt.imshow(img2, cmap='gray')  # Affichage

<matplotlib.image.AxesImage at 0x7f7c1655f040>

Conclusion de la partie

Vous savez maintenant:

• ce qu'est un signal et une fréquence d'échantillonnage
• lire, écrire, afficher des sons et des images en Python
• convertir des temps en indices entiers et vice-versa

On peut maintenant aborder des outils pour manipuler ces données

Introduction à la transformée de Fourier discrète

Transformée de Fourier discrète (DFT):

• notion fondamentale pour les signaux et pour l'UE SOAP
• outil pour comprendre les notions d'échantillonnage et de filtrage
• plusieurs aspects: algèbre linéaire, analyse spectrale

La DFT d'un vecteur $x$ de taille $N$ est un vecteur $X\in\mathbb{C}^N$ tel que $$ X = U^* x$$ où $U^*$ est une matrice $N\times N$ telle $U^*[k, n]= e^{-i2\pi\frac{kn}{N}}$.

C'est simplement un mapping linéaire!

Inversement, $x = \frac{1}{N} U X$.

De l'algèbre linéaire aux sinusoïdes...

$$ X = U^* x \text{ avec } U^*[k, n]= e^{-i2\pi\frac{kn}{N}}$$

$U^*=\begin{bmatrix}u_0^*\\ \vdots\\ u_{N-1}^*\end{bmatrix}$ est la transposée conjugée de $U = [u_0, \ldots, u_{N-1}]$ avec $u_k\in\mathbb{C}^N$ tel que $$u_k[n] = e^{i2\pi\frac{kn}{N}}=\cos(2\pi\frac{kn}{N}) + i \sin(2\pi\frac{kn}{N}).$$

• Les colonnes de $U$ forment une base orthogonale.
• $X$ contient les coefficients de la décomposition de $x$: $X[k]=<u_k, x>=u_k^* x$
• La DFT inverse donne la décomposition de $x$ dans cette base: $x = \sum_{k=0}^N X[k] u_k$

$\rightarrow$ $X[k]$ est le produit scalaire entre le signal et une sinusoïde complexe $u_k$!

Notion de fréquence

from numpy import pi

Signaux périodiques

Question: donnez la définition d'une fonction périodique

sound_filename = sound_path + '16.wav'
fs, x = wavfile.read(sound_filename)
plot_sound(x[:, 0], fs)
plt.xlim(6, 6.01)

(6.0, 6.01)

Quelles sont les notions et méthodes pour manipuler des signaux périodiques?

Sinusoïdes

La fonction sinus: $2\pi$-périodique, impaire.

x_range = np.linspace(-4*pi, 4*pi, 1000)
plt.plot(x_range, np.sin(x_range))

x_range = np.linspace(0, 2*pi, 1000)
plt.plot(x_range, np.sin(x_range))

plt.ylabel('sin(x)')
plt.xlabel('x')

Text(0.5, 0, 'x')

Remarque: la fonction cosinus est une translation de $\frac{\pi}{2}$ de sinus $$\forall x, \cos\left(x\right) = \sin\left(x+\frac{\pi}{2}\right)$$

# Vérification: superposition des deux fonctions
x_range = np.linspace(-5*pi, 5*pi, 1000)
plt.plot(x_range, np.cos(x_range), label='$\cos(x)$')
plt.plot(x_range, np.sin(x_range+pi/2), linestyle='-.', label='$\sin(x+\pi/2)$')
plt.xlabel('x')
plt.legend()

<matplotlib.legend.Legend at 0x7f7c04776640>

Une sinusoide de fréquence $f_0$ (en Hz) et de phase initiale $\varphi_0\in\mathbb{R}$ est définie pour tout temps $t$ (en secondes) par : $$ t \mapsto \sin(2\pi f_0 t + \varphi_0) $$ La sinusoide est de période $T_0=\frac{1}{f_0}$.

fs = 8000
duration = 2
t_range = np.arange(0, duration, step=1/fs)
f0 = 0.5
phi0 = 0
x0 = np.sin(2 * pi * f0 * t_range + phi0)

print(x0.shape)
plot_sound(x0, fs)

(16000,)

En faisant varier les fréquences:

for f in [0.25, 1, 4]:
    x = np.sin(2 * pi * f * t_range)
    plot_sound(x, fs, label=f'{f} Hz')
plt.legend()

<matplotlib.legend.Legend at 0x7f7c050abf70>

En faisant varier les phases initiales:

f0 = 0.5
for d in [4, 2, 1, 0.5]:
    x = np.sin(2 * pi * f0 * t_range + pi/d)
    plot_sound(x, fs, label=f'$\phi_0=\pi/{d}$')
plt.legend()

<matplotlib.legend.Legend at 0x7f7c050ab4c0>

Écoutons!

fs = 8000
for f in [100, 440, 1000, 3000]:
    x = np.sin(2 * pi * f * t_range)
    wavfile.write(filename=f'sin{f}.wav', data=x, rate=fs)

Les sinusoïdes complexes

Cas continu: sinusoïde complexe de fréquence $f_0 \in \mathbb{R}$ et phase initiale $\varphi_0 \in \left[0, 2\pi\right]$ $$ x(t) = e^{i (2\pi f_0 t + \varphi_0)} = \cos\left(2\pi f_0 t + \varphi_0\right) + i \sin\left(2\pi f_0 t + \varphi_0\right) = \left(e^{i 2\pi f_0}\right)^t e^{i \varphi_0}$$

Cas discret: on considère les instants discrets $t_n = \frac{n}{f_s}$ $$ x[n] = e^{i (2\pi f_0 \frac{n}{f_s} + \varphi_0)} = \cos\left(2\pi f_0 \frac{n}{f_s} + \varphi_0\right) + i \sin\left(2\pi f_0 \frac{n}{f_s} + \varphi_0\right) =\left(e^{i 2\pi \frac{f_0}{f_s}}\right)^n e^{i \varphi_0} $$

fs = 1000
f0 = 50
n_range = np.arange(1000)
z = np.exp(1j * 2 * pi * f0 * n_range / fs)
plt.figure(figsize=(5,5))
plt.plot(np.real(z), np.imag(z), '.')
plt.xlabel('Partie réelle')
plt.ylabel('Partie imaginaire')

plt.figure()
plt.plot(n_range/fs, np.real(z), '.-', label='Partie réelle')
plt.plot(n_range/fs, np.imag(z), '.-', label='Partie imaginaire')
plt.legend()
plt.xlim([0, 1/f0])
plt.xlabel('$t$')

plt.figure()
plt.plot(n_range/fs, np.real(z), label='Partie réelle')
plt.plot(n_range/fs, np.imag(z), label='Partie imaginaire')
plt.legend()
plt.xlabel('$t$')

Text(0.5, 0, '$t$')

Notion de fréquence normalisée/réduite

Conversion d'une fréquence $f$ en fréquence normalisée ou réduite $\nu \in \left[0, 1\right]$: $$\nu =\frac{f}{f_s}$$

Dans le cas de la sinusoïde discrète précédente, on a $\nu_0 n = f_0 \frac{n}{f_s} = f_0 t_n$ et on obtient $$ x[n] = e^{i (2\pi \nu_0 n + \varphi_0)} = \cos\left(2\pi \nu_0 n + \varphi_0\right) + i \sin\left(2\pi \nu_0 n + \varphi_0\right) =\left(e^{i 2\pi \nu_0}\right)^n e^{i \varphi_0}$$

Avantage: on ne considère que l'indice $n$ et une fréquence entre 0 et 1 (on ne considère plus le temps ni la fréquence d'échantillonnage, ni leurs unités s et Hz). Utile quand on veut simplifier certains calculs.

Exemple: construire un signal de taille 100 contenant une sinusoide avec 5 périodes

x_len = 100
nu = 5 / x_len  # inverse de la periode x_len / 5 = 100 / 5 = 20
x = np.sin(2 * pi * nu * np.arange(x_len))
plt.plot(x)
plt.title(rf'Sinusoïde de fréquence normalisée $\nu={nu}$')

Text(0.5, 1.0, 'Sinusoïde de fréquence normalisée $\\nu=0.05$')

Exercice: on considère la sinusoïde complexe discrète $x$ de fréquence normalisée $\nu_0 = 1/8$ et de phase initiale $\varphi_0$.

• Dans le plan complexe,
  • dessinez le cercle trigonométrique
  • placez $x[0]$ et faites apparaitre $\varphi_0$
  • placez $x[1], x[2], \ldots$
  • combien y a-t-il de points distincts et pourquoi?
• Sur un autre graphique, tracez la partie réelle et la partie imaginaire de $x$.

Conclusion de la partie

Vous avez maintenant une connaissance de l'élément de base pour l'analyse de signaux, la sinusoïde. Vous savez notamment:

• définir mathématiquement une sinusoïde réelle ou complexe
• interpréter les notions de fréquence et de phase initiale d'une sinusoïde
• représenter graphiquement ces signaux

On peut maintenant faire de l'analyse spectrale...

Analyse spectrale: la transformée de Fourier discrète

Représentation fréquentielle: intuition, exemples

Représentation fréquentielle d'une sinusoide

from soap_utils import plot_spectrum
fs = 8000
f_list = [200, 1000]
for f in f_list:
    x = np.sin(2 * pi * f * t_range)

    plt.figure()
    plt.subplot(121)
    plot_sound(x, fs, label=f'{f} Hz')
    plt.xlim(0, 1 / np.min(f_list))
    plt.legend()

    plt.subplot(122)
    plot_spectrum(x, fs, label=f'{f} Hz')
    plt.ylabel('Spectre')

Représentation fréquentielle d'un son

Dans un son, cherchons une note puis affichons un extrait et sa représentation fréquentielle

sound_path = '../../data/sons_int/'
i_sound = 35
sound_filename = sound_path + f'{i_sound:02d}.wav'
fs, x = wavfile.read(sound_filename)
x = x[:, 0]
y = x[int(0.6 * fs):int(0.7 * fs)]

plt.figure()
plt.subplot(121)
plot_sound(x, fs)
plt.title(f'Son {i_sound} entier')
plt.subplot(122)
plot_sound(x, fs)
plt.xlim(0.5, 2)
plt.title('Zoom sur un note')

plt.figure()
plt.subplot(121)
plot_sound(y, fs)
plt.title(f'Extrait de {i_sound} (100 ms)')
plt.subplot(122)
plot_spectrum(y, fs)
plt.title("Spectre de l'extrait")

Text(0.5, 1.0, "Spectre de l'extrait")

Les signaux périodiques et leur spectre harmonique

Exercice (TP): analyser le spectre de sons réels.

1. Pour les sons 08 à 14, faites comme précédemment en localisant une note et en analysant un extrait de 100 ms. Vous pouvez localiser l'extrait avec Audacity et ne pas afficher le son entier pour économiser de la mémoire : affichez uniquement l'extrait et de son spectre (vous pouvez faire un zoom sur les fréquences comprises entre $0$ et $f_s/2$ car les spectres sont symétriques).
2. Ces notes ont une hauteur précise : on parle de fréquence fondamentale $f_0$ de la note. La fréquence fondamentale est l'inverse de la période, que vous pouvez mesurer sur l'extrait affiché. Sur le spectre, les pics sont localisés aux multiples $f_n=n f_0$ de $f_0$, on peut donc aussi mesurer $f_0$ sur le spectre. Vous pouvez mesurer ces hauteurs. Pour convertir une fréquence fondamentale en nom de note (en notation anglaise: A pour la, B pour si, C pour do, etc., avec le numéro de l'octave), utilisez la fonction librosa.core.hz_to_note du module librosa; par exemple, $440$ Hz correspond au La 4 (A4).

# TODO

Transformée de Fourier discrète (DFT)

Intro: différentes transformées pour différents cas

Plusieurs variantes selon le type de signal: discret/continu, périodique ou pas.

Série de Fourier:

• pour des fonctions $T$-périodiques $x(t), t\in\mathbb{R}$
• définition: $c_k = \frac{1}{T}\int_{0}^T x(t) e^{-i 2\pi \frac{k}{T} t} dt$ pour $k\in\mathbb{Z}$
• fréquences $f_k=\frac{k}{T}$, multiples de $f_0=\frac{1}{T}$
• inverse: $x(t) = \sum_{k=-\infty}^{+\infty} c_k e^{i 2\pi \frac{k}{T} t}$

Transformée de Fourier

• pour des fonctions intégrables $x(t), t\in\mathbb{R}$
• définition $X(f) = \int_{-\infty}^{+\infty} x(t) e^{-i 2\pi f t} dt$ pour $f\in\mathbb{R}$
• fréquences $f\in\mathbb{R}$
• inverse $x(t) = \int_{-\infty}^{+\infty} X(f) e^{i 2\pi f t} dt$

Transformée de Fourier Discrète:

• pour des signaux $x[n], n\in\lbrace 0, \ldots, N-1 \rbrace$ de taille $N$
• définition: $X[k] = \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{kn}{N}}$ pour $k\in\lbrace 0, \ldots, N-1 \rbrace$
• fréquences $f_k=\frac{k}{N}f_s$, discrétisation de $[0, f_s[$
• inverse: $x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k] e^{i 2\pi \frac{kn}{N}}$

Définition et propriétés de la DFT

Définition. Pour un signal $x\in\mathbb{C}^N$, la DFT $X\in\mathbb{C}^N$ est définie pour tout $k\in\lbrace 0, \ldots, N-1\rbrace$ par $$X[k] = \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{kn}{N}}$$ On a bien $X[k]=\left\langle u_k, x \right\rangle$ où $u_k\in\mathbb{C}^N$ est défini pour $n\in\lbrace 0, \ldots, N-1\rbrace$ par $u_k[n]=e^{i 2\pi \frac{kn}{N}}$, le produit scalaire de $u, v\in\mathbb{C}^N$ étant défini par $\left\langle u, v\right \rangle=\sum_{n=0}^{N-1}\overline{u}[n] v[n]=u^* v$, où $u^*$ est le vecteur transposé conjugué de $u$.

On note aussi que $$X[k] = \left\langle c_k, x \right\rangle - i \left\langle s_k, x \right\rangle$$ où $c_k[n]= \cos\left(2\pi \frac{kn}{N}\right)$ et $s_k[n]= \sin\left(2\pi \frac{kn}{N}\right)$.

Propriétés (à prouver en TD):

• La famille $\lbrace u_0, \ldots, u_{N-1}\rbrace$ est une base orthogonale de $\mathbb{C}^N$. Comment en faire une base orthonormée?
• Transformée inverse: $\forall n, x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i 2\pi \frac{kn}{N}}$.
• Périodicité de la DFT: si l'on étend la définition $X[k] = \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{kn}{N}} dt$ à $k\in\mathbb{Z}$, alors $X$ est $N$-périodique.
• Symétrie hermitienne dans le cas d'un signal réel: si $x\in\mathbb{R}^N$, alors $\forall k, X[N-k] = \overline{X}[k]$.

Fréquences négatives: en utilisant les propriétés de symétrie et de périodicité, on peut donc représenter le spectre centré en $k=0$, c'est-à-dire pour $k\in \left\lbrace \left\lfloor -N/2+1\right\rfloor, \ldots, 0, \ldots, \left\lfloor N/2\right\rfloor\right\rbrace$. Les fréquences négatives ($k<0$) correspondent aux fréquences initialement définies pour $k\in ]N/2, N[$.

Représentation graphique d'un spectre

On représente en général le module de la DFT, en décibels (dB), c'est-à-dire $$20 \log_{10} \left|X[k]\right|$$

Exercice: DFT d'une sinusoïde et d'un dirac

• Calculez la TFD du signal $x[n] = e^{i2\pi \frac{k_0 n}{N}}$ où $k_0$ est un entier.
• En utilisant le résultat précédent, calculez la TFD de $x[n] = \cos\left (2\pi \frac{k_0 n}{N}\right )$. Où sont les pics de cette TFD? Commentez le comportement des parties réelles et imaginaires de cette TFD.
• Calculez la TFD du signal $\delta_{n_0}[n] = \begin{cases}1 \text{ si } n=n_0\\0 \text{ sinon }\\\end{cases}$ où $n_0$ est un entier. Tracez le signal et son spectre.

sound_path = '../../data/sons_int/'
i_sound = 35
sound_filename = sound_path + f'{i_sound:02d}.wav'
fs, x = wavfile.read(sound_filename)
x = x[:, 0]
y = x[int(0.6 * fs):int(0.7 * fs)]

plt.figure()
plt.subplot(121)
plot_sound(y, fs)
plt.title(f'Extrait de {i_sound} (100 ms)')
plt.subplot(122)
plot_spectrum(y, fs)
plt.title("Spectre de l'extrait")

plt.figure()
plt.subplot(121)
plt.plot(np.abs(np.fft.fft(y)))
plt.title(f'Spectre : affichage linéaire')
plt.subplot(122)
plt.plot(20 * np.log10(np.abs(np.fft.fft(y))))
plt.title("Spectre en dB")

Text(0.5, 1.0, 'Spectre en dB')

Illustrations numériques des propriétés de la DFT

Base orthogonale

Vérifions que $U=\frac{1}{\sqrt{N}}[u_0, \ldots, u_{N-1}]$ est la matrice d'une base orthonormée.

n = 128  # dimension du signal
n_range = np.arange(n)
U = np.exp(1j * 2 * pi *np.outer(n_range, n_range / n)) / np.sqrt(n)

UU = np.conj(U).T @ U

plt.figure(figsize=(20, 8))
plt.subplot(221)
plt.imshow(np.real(UU))
plt.colorbar()
plt.title('$Re(U^*U)$')
plt.subplot(222)
plt.imshow(np.real(UU)-np.eye(n))
plt.colorbar()
plt.title('$Re(U^*U)-I_N$')
plt.subplot(223)
plt.imshow(np.imag(UU))
plt.colorbar()
plt.title('$Im(U^*U)$')
plt.subplot(224)
plt.plot(np.real(np.diag(UU)))
plt.title('$Re(diag(U^*U))$')

Text(0.5, 1.0, '$Re(diag(U^*U))$')

Affichons quelques vecteurs de la base:

for k in [0, 1, 2, 4, n//2-1]:
    plt.subplot(121)
    plt.plot(np.real(U[:, k]))
    plt.title('Partie réelle')
    plt.subplot(122)
    plt.plot(np.imag(U[:, k]), label=str(k))
    plt.title('Partie imaginaire')
plt.legend()

<matplotlib.legend.Legend at 0x7f7bfdcc2af0>

DFT d'un son:

sound_path = '../../data/sons_int/'
i_sound = 35
sound_filename = sound_path + f'{i_sound:02d}.wav'
fs, x = wavfile.read(sound_filename)
x = x[int(0.6 * fs):int(0.6 * fs)+n, 0]
X = np.conj(U).T @ x

plt.subplot(121)
plot_spectrum(x, fs)
plt.title('Calcul avec la fonction fft')
plt.subplot(122)
plt.plot(20 * np.log10(np.abs(X) * np.sqrt(n)))
plt.title('Calcul notre matrice U')

Text(0.5, 1.0, 'Calcul notre matrice U')

Périodicité et symétrie

Que se passe-t-il si l'on manipule des fréquences en dehors de $[0, f_s/2]$?

Les exemples suivants montrent que ces signaux numériques coïncident avec des signaux de fréquences dans $[0, f_s/2]$.

Périodicité du spectre

La DFT d'un signal de taille $N$ est $N$-périodique.

Conséquence: les versions numériques d'une sinusoïde à la fréquence $f_0$ et d'une autre à la fréquence $f_s+f_0$ sont égales. $$ \forall f_0, \forall n\in\mathbb{Z}, \sin\left(2\pi \left(f_s+f_0\right)t_n\right) = \sin\left(2\pi f_0 t_n\right) \text{ avec } t_n=\frac{n}{f_s}$$

fs = 8000

plt.figure(figsize=(20, 10))
f_list = [200, 1000]
for i in range(len(f_list)):
    fx = f_list[i]
    fy = fs + fx

    x = np.sin(2 * pi * fx * t_range)
    X = np.fft.fft(x)
    y = np.sin(2 * pi * fy * t_range)
    Y = np.fft.fft(y)

    print('erreur:', np.max(np.abs(x-y)))

    plt.subplot(len(f_list), 2, 2 * i + 1)
    plot_sound(x, fs, label=f'{fx} Hz')
    plot_sound(y, fs, label=f'{fy} Hz', linestyle='-.')
    plt.xlim(0, 1 / np.min(f_list))
    plt.legend()

    plt.subplot(len(f_list), 2, 2 * i + 2)
    plot_spectrum(x, fs, label=f'{fx} Hz')
    plot_spectrum(y, fs, label=f'{fy} Hz', linestyle='-.')
    plt.ylabel('Spectre')

erreur: 1.7631859650566897e-11
erreur: 1.762668600371223e-11

Symétrie hermitienne

$x\in\mathbb{R}^N$ est un signal réel ssi $ \forall k, X[k] = X[N-k]^*$

Conséquence: les versions numériques d'une sinusoïde à la fréquence $f_0$ et d'une autre à la fréquence $f_s-f_0$ sont égales, au signe près. $$ \forall f_0, \forall n\in\mathbb{Z}, \sin\left(2\pi \left(f_s-f_0\right)t_n\right) = -\sin\left(2\pi f_0 t_n\right) \text{ avec } t_n=\frac{n}{f_s}$$

fs = 8000

plt.figure(figsize=(20, 10))
f_list = [200, 1000]
for i in range(len(f_list)):
    fx = f_list[i]
    fy = fs - fx

    x = np.sin(2 * pi * fx * t_range)
    X = np.fft.fft(x)
    y = np.sin(2 * pi * fy * t_range)
    Y = np.fft.fft(y)
    
    print(np.max(np.abs(x+y)))
    
    plt.subplot(len(f_list), 2, 2 * i + 1)
    plot_sound(x, fs, label=f'{fx} Hz')
    plot_sound(y, fs, label=f'{fy} Hz', linestyle='-.')
    plt.xlim(0, 1 / np.min(f_list))
    plt.legend()

    plt.subplot(len(f_list), 2, 2 * i + 2)
    plot_spectrum(x, fs, label=f'{fx} Hz')
    plot_spectrum(y, fs, label=f'{fy} Hz', linestyle='-.')
    plt.ylabel('Spectre')

1.5680845510956942e-11
1.821205044844679e-11

Représentation alternative avec fréquences négatives

fs = 8000

plt.figure(figsize=(20, 10))
f_list = [200, 1000]
for i in range(len(f_list)):
    fx = f_list[i]
    fy = fs + fx

    x = np.sin(2 * pi * fx * t_range)
    X = np.fft.fft(x)

    plt.subplot(len(f_list), 3, 3 * i + 1)
    plot_sound(x, fs, label=f'{fx} Hz')
    plt.xlim(0, 1 / np.min(f_list))
    plt.legend()

    plt.subplot(len(f_list), 3, 3 * i + 2)
    plot_spectrum(x, fs, label='f{fx} Hz')
    plt.ylabel('Spectre')
    
    plt.subplot(len(f_list), 3, 3 * i + 3)
    plot_spectrum(x, fs, fft_shift=True, label=f'{fx} Hz')
    plt.ylabel('Spectre')

Algorithme rapide: fast Fourier transform (FFT)

Algorithme naïf $X=U \times x$ de complexité $\mathcal{O}(N^2)$ en temps et en espace.

Cooley–Tukey FFT (1965): algorithme rapide de complexité $\mathcal{O}(N\log(N))$ en temps et en espace.

Conclusion de la partie

Vous avez maintenant une connaissance approfondie de la transformée de Fourier discrète, qui permet d'analyser le contenu fréquentiel des sons et autres signaux 1D. Vous savez

• définir mathématiquement la DFT et interpréter cette définition en termes de produit scalaire avec des sinusoïdes
• un certain nombre de propriétés de la DFT, leur interprétation
• représenter graphiquement le spectre de signaux

On peut maintenant faire la même chose en 2D avec les images...

Fréquences spatiales et DFT pour les images

La DFT s'étend aux images:

• dimension des données: remplacer le temps par l'espace $\rightarrow$ notion de fréquence spatiale
• passer de signaux 1D à 2D $\rightarrow$ DFT 2D

DFT 2D d'une image

Définition. Soit $I \in \mathbb{R}^{N \times M}$ une image. La DFT 2D de $I$ est une matrice $S\in \mathbb{C}^{N \times M}$ obtenue en faisant la DFT de chaque ligne suivie de la DFT de chaque colonne de la matrice obtenue. Pour $k\in\lbrace 0, \ldots, N-1\rbrace$ et $l\in\lbrace 0, \ldots, M-1\rbrace$, on a $$S[k, l]=\sum_n \sum_m I[n,m] e^{-i2\pi\left(\frac{kn}{N}+\frac{lm}{M}\right)}$$

$\frac{k}{N}$ est la $k$-ième fréquence verticale, $\frac{l}{M}$ est la $l$-ième fréquence horizontale.

Propriétés:

• les propriétés de la DFT 1D s'étendent au cas 2D
• la DFT 2D peut aussi être obtenue en faisant la DFT des colonnes suivie de la DFT des lignes de la matrice obtenue.

Spectre 2D d'une image

from soap_utils import show_spectrum_2d

# Chargement d'une image
image_path = '../../data/images/'
image_filename = image_path + 'barbara.png'
img = plt.imread(image_filename)
plt.figure(figsize=(20, 10))
plt.imshow(img, cmap='gray')
plt.colorbar()
print(f'Image chargée: type {type(img)}, format {img.shape}, dtype {img.dtype}')
plt.figure(figsize=(20, 10))
show_spectrum_2d(img, db_scale=True)

Image chargée: type <class 'numpy.ndarray'>, format (512, 512), dtype float32

Spectre 2D d'un extrait d'image

# Extraction d'une partie
imgz = img[230:260, 290:350]

plt.figure(figsize=(15, 6))
plt.subplot(121)
plt.imshow(imgz, cmap='gray')
plt.title('Extrait')
plt.subplot(122)
show_spectrum_2d(imgz, db_scale=True)
plt.title('Spectre 2D')

Text(0.5, 1.0, 'Spectre 2D')

Les fréquences horizontales

k = 0
for l in [0, 1, 2, 5, 15, 25]:
    S = np.zeros((50, 50))
    S[k, l] = 1
    img = np.real(np.fft.ifft2(S))
    plt.figure(figsize=(20, 4))
    plt.subplot(131)
    plt.imshow(img)
    plt.title(str((k, l)))
    plt.subplot(132)
    plt.plot(img[0, :])    
    plt.subplot(133)
    show_spectrum_2d(img)
    plt.title('Spectre')

Les fréquences verticales

l = 0
for k in [0, 1, 2, 5, 15, 25]:
    S = np.zeros((50, 50))
    S[k, l] = 1
    img = np.real(np.fft.ifft2(S))
    plt.figure(figsize=(20, 4))
    plt.subplot(131)
    plt.imshow(img)
    plt.title(str((k, l)))
    plt.subplot(132)
    plt.plot(img[:, 0])    
    plt.subplot(133)
    show_spectrum_2d(img)
    plt.title('Spectre')

Les fréquences obliques

Un couple (k, l) d'indices correspondant à une fréquence verticale $\nu_v=\frac{k}{N}$ et une fréquence horizontale $\nu_h=\frac{l}{M}$ correspond à une fréquence oblique $\nu=\sqrt{\left(\nu_h^2+\nu_v^2\right)}$.

for k, l in [(1, 1), (1, 5), (-1, 5), (1, 15), (10, 15), (25, 25), (45, 1)]:
    S = np.zeros((50, 50))
    S[k, l] = 1
    img = np.real(np.fft.ifft2(S))
    plt.figure(figsize=(20, 4))
    plt.subplot(131)
    plt.imshow(img)
    plt.title(str((k, l)))
    plt.subplot(132)
    plt.plot(img[3, :])
    plt.plot(img[:, 4]) 
    plt.subplot(133)
    show_spectrum_2d(img)
    plt.title('Spectre')

Exercice: supprimer les rayures sur l'image barbara.png

On procèdera ainsi:

• chargez et affichez l'image $N\times M$ barbara.png
• on veut repérer la période des rayures à enlever
  • extrayez une zone avec des rayures, comme dans l'exemple déjà vu
  • affichez le spectre 2D de cette zone et repérez les fréquences horizontales $\nu_h$ et verticales $\nu_v$ des rayures. On rappelle qu'une fréquence d'indice $k$ a pour valeur $\frac{k}{N}$ où $N$ est le nombre de pixels dans la direction en question.
  • calculez ensuite la fréquence oblique correspondante $\nu=\sqrt{\nu_h^2+\nu_v^2}$
• calculez le spectre 2D de l'image entière, puis
  • annulez les colonnes correspondant à toutes les fréquences horizontales supérieures à $\nu$, et reconstruisez l'image (np.fft.ifft2 suivi de np.real) et affichez-là: que constatez-vous?
  • recommencez l'étape précédente en annulant les lignes correspondant aux fréquences verticales supérieures à $\nu$
  • recommancez en annulant à la fois les lignes et les colonnes repérées
  • que faudrait-il vraiment faire pour annuler toutes les composantes de fréquence supérieure à $\nu$?

Remarque: n'hésitez pas à prendre des valeurs assez faibles pour $\nu_h$ et $\nu_v$ (le début du pic plutôt que le sommet) pour enlever ensuite toute la partie correspondante.

# TODO

Autre application: échantillonnage k-space en imagerie

Voir k-Space tutorial: an MRI educational tool for a better understanding of k-space

Conclusion de la partie

Vous maitrisez maintenant la version 2D de la DFT:

• sa définition
• les notions de fréquences verticales, horizontales et obliques
• l'interprétation du spectre 2D d'images

On peut maintenant passer à un résultat majeur: le théorème d'échantillonnage de Shannon-Nyquist.

Comment échantillonner un signal ou une image?

Les fréquences discrètes dans la DFT

Avec une fréquence d'échantillonnage $f_s$, on a $$u_k[n] = e^{i2\pi\frac{kn}{N}} = e^{i2\pi f_k t_n} $$ avec

• $f_k=\frac{k}{N}f_s$ en Hz
• $t_n=\frac{n}{f_s}$ en secondes

$\rightarrow$ Les fréquences positives de la DFT d'un signal réel sont comprises entre $f_0=0$ et $f_{N/2}=\frac{f_s}{2}$.

Théorème d'échantillonnage de Shannon-Nyquist

Énoncé. Tout signal peut être échantillonné régulièrement à une fréquence d'échantillonnage $f_s$ sans perte d'information si son contenu fréquentiel est compris entre $0$ et $f_s/2$.

Autrement dit, la fréquence d'échantillonnage $f_s$ doit être au moins deux fois plus élevée que la fréquence maximale présente dans un signal pour ne pas perdre d'information ni déteriorer le signal.

Utilités:

• numériser un signal continu
• sous-échantillonner un signal numérique

En pratique, on applique d'abord un filtrage passe-bas qui élimine tout le contenu fréquentiel au-delà de $f_s/2$ avant l'échantillonnage.

Si l'on ne respecte pas les conditions du théorème, on retrouve entre $0$ et $f_s/2$ l'énergie des fréquences qui étaient en dehors de cet intervalle: c'est le phénomène de repliement spectral.

Question: sachant que l'oreille humaine perçoit des fréquences comprises entre $20$ et $20 000$ Hz, expliquez la fréquence d'échantillonnage des CD.

Formule de reconstruction

Sous les hypothèses du théorème, pour reconstruire le signal $x(t)$ en fonction de son échantillonnage $x[n]$, on a $$\forall t \in \mathbb{R}, x(t) = \sum_{n} x[n] \times sinc\left(f_s t - n\right)$$ où $sinc(t)=\frac{\sin\left(\pi t\right)}{\pi t}$ pour $t\neq 0$ et $sinc(0)=1$.

Comment sous-échantillonner

Ce qu'il ne faut pas faire...

Plusieurs exemples où l'on sous-échantillonne avec du repliement spectral...

# Exemple de repliement

fs = 8000
duration = 1
t_range = np.arange(0, duration, step=1/fs)

for f0 in [100, 2100]:
    x0 = np.sin(2 * pi * f0 * t_range)

    downsampling_factor = 4
    fs_down = fs / downsampling_factor
    x0_down = x0[::downsampling_factor]

    print(x0.shape)

    plt.figure(figsize=(20, 4))
    plt.subplot(121)
    plot_sound(x0, fs, label=f'son original ($f_s={fs}$)', marker='.')
    plot_sound(x0_down, fs_down, label=f'son sous échantillonné ($f_s={fs_down}$)', marker='.', linestyle='-.')
    plt.xlim(0, 5/f0)
    plt.legend()
    plt.title(f'$f_0={f0}$')
    plt.subplot(122)
    plot_spectrum(x0, fs, fft_shift=True, label='son original')
    plot_spectrum(x0_down, fs_down, fft_shift=True, label='son sous échantillonné', linestyle='-.')

(8000,)
(8000,)

# Dans cet exemple, on "annule une fréquence f1" via une fréquence f2

fs = 8000
duration = 1
t_range = np.arange(0, duration, step=1/fs)

# Création d'un signal avec deux sinusoïdes
f1 = 100
f2 = 2100
x = np.sin(2 * pi * f1 * t_range) + np.sin(2 * pi * f2 * t_range + pi)

# Sous-échantillonnage
downsampling_factor = 4
fs_down = fs / downsampling_factor
x_down = x[::downsampling_factor]

plt.figure(figsize=(20, 4))
plt.subplot(121)
plot_sound(x, fs, label=f'son original ($f_s={fs}$)')
plt.xlim(0, 2/f1)
plt.legend()
plt.subplot(122)
plot_spectrum(x, fs, fft_shift=True, label='son original')

plt.figure(figsize=(20, 4))
plt.subplot(121)
plot_sound(x, fs, label=f'son original ($f_s={fs}$)')
plot_sound(x_down, fs_down, label=f'son sous échantillonné ($f_s={fs_down}$)')
plt.legend()
plt.subplot(122)
plot_spectrum(x, fs, fft_shift=True, label='son original')
plot_spectrum(x_down, fs_down, fft_shift=True, label='son sous échantillonné', linestyle='-.')

Ce qu'il faut faire:

from scipy.signal import resample

fs = 8000
duration = 1
t_range = np.arange(0, duration, step=1/fs)

for f0 in [100, 2100]:
    x0 = np.sin(2 * pi * f0 * t_range)

    downsampling_factor = 4
    fs_down = fs / 4
    x0_down = resample(x0, int(x0.size / downsampling_factor))

    print(x0.shape)

    plt.figure(figsize=(20, 4))
    plt.subplot(121)
    plot_sound(x0, fs, label=f'son original ($f_s={fs}$)', marker='.')
    plot_sound(x0_down, fs_down, label=f'son sous échantillonné ($f_s={fs_down}$)', marker='.', linestyle='-.')
    plt.xlim(0, 5/f0)
    plt.legend()
    plt.title(f'$f_0={f0}$')
    plt.subplot(122)
    plot_spectrum(x0, fs, fft_shift=True, label='son original')
    plot_spectrum(x0_down, fs_down, fft_shift=True, label='son sous échantillonné', linestyle='-.')

(8000,)
(8000,)

Exercice: repliement fréquentiel. Cet exercice peut être fait entièrement en Python ou à l'aide des logiciels d'édition Audacity et The Gimp (pour certains sous-échantillonnages)

• Affichez les DFT 1D des sons et les DFT 2D des images pour sélectionner les données contenant des fréquences élevées.
• Pour chaque fichier sélectionné, produisez des versions sous-échantillonnées avec et sans repliement spectral
• Affichez côte à côte (ou en superposition dans le cas des sons): le spectre du fichier original et les spectres des deux versions sous-échantillonnées pour visualiser le repliement
• Comparez les trois versions d'un point de vue perceptif : visuellement pour les images, auditivement pour les sons
• Question subsidiaire: pour exhiber le défaut, créez une nouvelle variante contenant uniquement la partie repliée du spectre puis visualisez l'image obtenu/écoutez le son obtenu.